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Abstract 

We apply the recently devised quasi-stationary simulation method to study the lifetime and 
order parameter of the contact process in the subcritical phase. This phase is not accessible to 
other methods because virtually all realizations of the process fall into the absorbing state before 
the quasi-stationary regime is attained. With relatively modest simulations, the method yields an 
estimate of the critical exponent vu with a precision of 0.5%. 
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I. INTRODUCTION 

Stochastic processes with an absorbing state arise frequently in statistical physics 
epidemiology |3|, and related fields. Phase transitions to an absorbing state in spatially 
extended systems, exemplified by the contact process 0,0], are currently of great interest 
in connection with self-organized criticality p, the transition to turbulence p|, and issues 



of universality in nonequilibrium critical phenomena 

Systems exhibiting a phase transition to absorbing an state possess (for appropriate values 
of the control parameter), an active (nonabsorbing) stationary state in the infinite-size limit. 
But if the system size is finite, the process must eventually end up in the absorbing state. 
The quasi- stationary (QS) distribution for such a system provides a wealth of information 
about its behavior. (Since the only true stationary state for a finite system is the absorbing 
one, simulations of "stationary" properties of models with an absorbing state in fact study 
the quasi-stationary regime.) 

Conventional Monte Carlo methods entail a somewhat complicated procedure for deter- 
mining QS properties: many independent realizations are performed (using the same initial 
configuration), and the mean (j)(t) of some property (for example the order parameter) is 
evaluated over the surviving realizations at time t. At short times times 4>(t) exhibits a 
transient as it relaxes toward the QS regime; at long times it fluctuates wildly as the sur- 
viving sample decays. In the supercritical phase (where an infinite system survives forever) 
one is normally able to identify an intermediate regime free of transients and with limited 
fluctuations, which can be used to estimate the QS value of d>. (Even in this case the method 
requires careful scrutiny of the data and is not always free of ambiguity [12|].) Deep in the 
subcritical phase, however, conventional simulations are impractical, as nearly all realiza- 
tions fall into the absorbing state before the quasi-stationary regime is attained. (In other 
words, one no longer has a separation of the time scales for relaxation to the QS state and 
for its survival. The effective potential governing the process no longer possesses a local 
minimum away from the absorbing state.) 

We recently devised a simulation method that yields quasi-stationary properties directly 
jl^l ]. We showed that the method reproduces known scaling properties of the contact process 
(CP); in [1J| it was used to study the static correlation function of the model. In studies of 
the critical point, the QS simulation method requires an order of magnitude less computer 
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time than conventional simulations, to obtain results of comparable precision. Here we use 
the method to study the lifetime and the order parameter of the one-dimensional contact 
process in the subcritical phase. 

In the following section we review the method and define the contact process. Then in 
Sec. Ill we present our results for the lifetime of the CP on a ring, in the subcritical phase. 
We summarize our findings in Sec. IV. 



II. BACKGROUND 



We begin by reviewing the definition of the quasi-stationary distribution. Consider a 
continuous-time Markov process X t taking values n = 0, 1, 2, S, with the state n = 
absorbing. (In any realization of the process, if X t = 0, then X t > must be zero for all 
subsequent times, t' > t.) The transition rates w m ^ n (from state n to state m) are such that 
w m ,o = 0, Vm > 0. Let p n (t) denote the probability of state n at time t, given some initial 
state Xq. The survival probability P s (t) = Y2 n >iP n (t) * s ^ ne probability that the process 
has not become trapped in the absorbing state up to time t. We suppose that as t — > oo 
the p n , normalized by the survival probability P s (t), attain a time-independent form. The 
quasi- stationary distribution p n is then defined via 

t->oo f s (t) 

with p = 0. The QS distribution is normalized so: 

Z>» = 1 - (2) 

n>l 

n 

The basis for our simulation method is the observation pj| that the QS distribution 
represents a stationary solution to the equation, 

^ = -w n q n + r n + r q n (n > 0) , (3) 

where w n = Ylim w m,n is the total rate of transitions out of state n, and r n = J2 m w n,m ( lm 
is the flux of probability into this state. To see this, consider the master equation (Eq. 
without the final term) in the QS regime. Substituting q n (t) = P s (t)p n , and noting that 
in the QS regime dP s /dt = —r = —P s J2m w o,™p m i we see that the r.h.s. of Eq. (jSJ) is 
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identically zero if q n = p n for n > 1. The final term in Eq. (jSJ) represents a redistribution of 
the probability tq (transfered to the absorbing state in the original master equation), to the 
nonabsorbing subspace. Each nonabsorbing state receives a share equal to its probability. 

Although Eq. (J3J) is not a master equation (it is nonlinear in the probabilities q n ), it 
does suggest a simulation scheme for sampling the QS distribution. In a Monte Carlo 
simulation one generates a set of realizations of a stochastic process. In what follows we 
call a simulation of the original process X t (possessing an absorbing state) a conventional 
simulation. We define a related process, X*, whose stationary probability distribution is 
the quasi-stationary distribution of X t . (Note that in order to have a nontrivial stationary 
distribution, X t * cannot possess an absorbing state.) The probability distribution of X* is 
governed by Eq. 0, which implies that for n > (i.e., away from the absorbing state), the 
evolution of X£ is identical to that of X t . (Since Eq. (jSJ) holds for n > 0, the process X 4 * 
must begin in a non-absorbing state.) When X t enters the absorbing state, however, X 4 * 
instead jumps to a nonabsorbing one, and then resumes its "usual" evolution (i.e., with the 
same transition probabilities as X t ), until such time as another visit to the absorbing state 
is imminent. The probability that X 4 * jumps to state j (when a visit to state is imminent), 
is the QS probability py 

A subtlety associated with this procedure is that the QS distribution is needed to deter- 
mine the evolution of A 4 * when X t visits the absorbing state. Although one has no prior 
knowledge of the QS distribution p n , one can, in a simulation, use the history X* (0 < s < t) 
up to time t, to estimate the p n . (There is good evidence, after all, that the surviving sample 
in a conventional simulation converges to the QS state after an initial transient.) In practice 
this is accomplished by saving (and periodically updating) a sample ni, n 2 , Um of the 
states visited. As the evolution progresses, X* will visit states according to the QS distribu- 
tion. We therefore update the sample {ni, ri2, Um} by occasionally replacing one of these 
configurations with the current one. In this way the distribution for the process X* (and 
the sample drawn from it), will converge to the QS distribution (i.e., the stationary solution 
of Eq. (jnj)) at long times. Summarizing, the simulation process A t * has the same dynamics 
as X t , except that when a transition to the absorbing state is imminent, X* is placed in a 
nonabsorbing state, selected at random from a sample over the history of the realization. In 
effect, the nonlinear term in Eq. © is represented as a memory in the simulation. 
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A. The contact process 

To explain how our method works in practice, we detail its application to the contact 
process (CP) DDD . In the CP, each site % of a lattice is either occupied (<7j(t) = 1), or 
vacant (o~i(t) = 0). Transitions from <jj = 1 to Oi = occur at a rate of unity, independent 
of the neighboring sites. The reverse transition is only possible if at least one neighbor is 
occupied: the transition from o~.j = to ^ = 1 occurs at rate Ar, where r is the fraction of 
nearest neighbors of site i that are occupied; thus the state <jj = for all i is absorbing. (A 
is a control parameter governing the rate of spread of activity.) 

Although no exact results are available, the CP has been studied intensively via series 
expansion and Monte Carlo simulation. The model has attracted much interest as a proto- 
type of a nonequilibrium critical point, a simple representative of the directed percolation 
(DP) universality class. Since its scaling properties have been discussed extensively (si 
we review them only briefly. The best estimate for the critical point in one dimension is 
A c = 3.297848(20), as determined via series analysis As the critical point is approached, 
the correlation length £ and correlation time r diverge, following £ oc |A| _I/± and r oc |A| _!/ ii, 
where A = (A — A c )/A c is the relative distance from the critical point. The order parameter 
(the fraction of active sites), scales as p oc A 13 for A > 0. 

Two characteristic times, tq and tl, may be identified in the contact process. The first 
is a relaxation time that governs the decay of temporal correlations in the stationary state: 
C(t) = (p(to)p(to + t)) — p 2 ~ e~ i//Tc . The second is a lifetime, determining the asymptotic 
decay of the survival probability (starting from a spatially homogeneous initial condition) 
via P(t) ~ e~ l l TL . The two characteristic times exhibit the same scaling properties in the 
critical region. In the supercritical or active phase (A > 0), the lifetime grows exponentially 
with system size and A, while tq remains finite. Our interest here is in the lifetime (denoted 
simply as r) of the quasi-stationary state. In the QS probability distribution there is a 
nonzero flux of probability to the absorbing state, 



r = w tpi (4) 

where pi is the QS probability of the configuration with exactly one active site and u>oi = 1 is 
the transition rate from this configuration to the absorbing state; the QS lifetime r = 1/tq. 
(In QS simulations we take r to be the mean time between successive attempts to visit to 
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the absorbing state.) 

An important point in interpreting our simulation results concerns the finite-size scaling 
(FSS) behavior of r. According to the usual FSS hypothesis [l^l, finite-size corrections to 
critical properties are functions of the ratio or, equivalently, of the quantity AL 1 ^ U± . 
The lifetime is therefore expected to follow (in the subcritical phase, A < 0), 

r(A,L) = \A\- U ^(\A\L 1 /^) (5) 

where the scaling function J-(x) oc x v w for small x (so that r does not diverge in a finite 
system), while in the opposite limit T — ► J-"o, a constant. The scaling hypothesis leads to 
the familiar relation, r(0, L) ~ L v \\/ V± at the critical point, and suggests that we attempt to 
collapse data for diverse system sizes by plotting A"Ht versus A* = AL 1 ^ U± . For the order 
parameter the expected finite-size scaling form is , 

p(A,L) = \Af'Jl(L 1/ ^A) . (6) 

In the subcritical regime, the order parameter must decay to zero oc L^ 1 as L — > oo, for any 
A < 0, so that 1Z(x) ~ l^l - ^ 1 " as x —>■ — oo. On the other hand, for A = and L finite, p 
must be nonzero and nonsingular, implying TZ(x) ~ x~^ for x — > 0. Thus p ~ |A|^ _i/± for 
A large and negative. 



III. SIMULATION RESULTS 

In the QS simulations we use a list of size M = 2 xlO 3 - 10 4 , depending on the lattice size. 
The process is simulated in 15 realizations, each of 5x 10 8 time steps. As is usual, annihilation 
events are chosen with probability 1/(1 + A) and creation with probability A/(l + A). A site 
i is chosen from a list of currently occupied sites, and, in the case of annihilation, is vacated. 
In a creation event, a nearest-neighbor of site % is selected at random and, if it is currently 
vacant, it becomes occupied. The time increment associated with each event is At = l/N occ , 
where N occ is the number of occupied sites just prior to the attempted transition 0. 

In the initial phase of the evolution, the list of saved configurations is augmented whenever 
the time t increases by one, until a list of M configurations has been accumulated. From then 
on, we update the list (replacing a randomly selected entry with the current configuration), 
with a certain probability p rep , whenever t advances by one unit. A given configuration 
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therefore remains on the list for a mean time of M/p rep . (Values of p rep in the range 
10~ 3 — 10~ 4 are used.) 

Fig. 1 shows the QS lifetime r as a function of A, for lattice sizes L = 20, 40, 80,. ..,2560. 
For the larger system sizes, power-law dependence on A is evident, before the lifetime 
saturates (at very small values of |A|), due, as anticipated, to finite-size effects. In Fig. 2 
these data are collapsed using the known values of the DP exponents ^j], v± = 1.09684(6) 
and v\\ = 1.73383(3). The data collapse is quite good for the larger system sizes. A least- 
squares fit to the data in the linear portion of the graph (|A| > 0.02) yields a slope of 
-1.738(12), in reasonable agreement with the accepted value for u\\. An alternative method 
for analyzing the data is to estimate, for each value of A, the infinite-size limiting value of 
the lifetime by plotting r(A, L) versus 1/L (see Fig. 3). Plotting the resulting estimates on 
log scales, we find a slope of -1.735(9), again in good agreement with the standard value for 
the exponent v\\. 

We verified that the distribution px of the lifetime (that is, of the time interval between 
successive attempts to visit the absorbing state) is exponential: 

p-t/r 

Mt) = (7) 

T 

This is as expected: in the quasi-stationary state the transition rate to the absorbing state 
is time-independent. 

Turning to the order parameter, we see from Eq. (JHJ) that a plot of p* = |A| _/3 p versus A* 
should exhibit a data collapse. The asymptotic behavior of the scaling function 1Z implies 
that p* oc (A*)~ u± for large A*. These scaling properties are verified in Fig. 4. (We use the 
accepted value f3 = 0.2765. While scaling of the order parameter in the subcritical regime 
was verified in Ref here we are able to extend the range of A* by an order of magnitude, 
using the QS simulation technique.) 

In Fig. 5 we show the probability distribution p(p) of the order parameter in a large 
system (L = 2560) for three values of A. As expected, the distribution broadens and shifts 
toward larger values of p as A approaches the critical value. The probability distribution 
follows, to a good approximation, the scaling form 

p(p) = yfipHp)) (8) 
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where (p) = p(A, L) is the mean value and V is a scaling function. This is verified (Fig. 
6) by plotting p* = (p)p(p) versus p/(p). We see that the scaling function P attains its 
maximum near p* = 0.6, and that it falls off rapidly as p* — > 0. On the other side of the 
maximum it exhibits a roughly exponential tail. 

IV. SUMMARY 

We use the quasi-stationary simulation method to study the lifetime and order parameter 
of the one-dimensional contact process in the subcritical phase. Our results confirm the 
expected scaling properties of the lifetime of the quasi-stationary state, and of the order 
parameter. The QS simulation method is the first to allow such an analysis deep in the 
subcritical regime. With a rather modest expenditure of computer time, the approach 
yields an estimate of the critical exponent uu that agrees with the accepted value to within 
uncertainty, with a precision of about 0.5%. Analysis of quasi-stationary properties in the 
subcritical regime should therefore be a useful tool in the study of absorbing-state phase 
transitions. 
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FIGURE CAPTIONS 

FIG. I. QS survival time r as a function of A, for lattice sizes L = 20, 40, 80,.. .,2560. The 
slope of the straight line is -1.74. 

FIG. 2. Data of Fig.l plotted in terms of the scaling variables A* = L 1 ^ U± \A\ and r* = L~ z t, 
with z = v\\jv±_. The slope of the straight line is -1.734. 

FIG. 3. QS survival time r as a function of A, for |A| = 

0.3,0.2,0.1,0.05,0.02,0.01,0.005,0.002,0.001,0.0005, from bottom to top (lattice sizes 
L = 20, 40, 80,.. .,2560). Lifetime by plotting r(A, L) versus l/L. Inset: infinite-size 
limiting estimates of the lifetime as function of |A|, plotted on log scales. 
FIG. 4. Scaled QS order parameter p* = \A\~Pp versus A* = L l ' v ^\A\ for L = 320 (open 
squares), L = 1280 (diamonds) and L = 2560 (filled squares). The slope of the straight line 
is -1.1097. 

FIG. 5. QS probability distribution of the order parameter in the subcritical regime, for 
L = 2560 and (left to right) A/A c = 0.8, 0.9 and 0.95. 

FIG. 6. Scaling plot of the data shown in Fig. 5; A/A c = 0.8 (squares), 0.9 (circles) and 0.95 
(triangles). 
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